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Abstract. Defects are believed to play a fundamental role in the supersolid 
state of 4 He. We have studied solid 4 He in two dimensions (2D) as function 
of the number of vacancies n v , up to 30, inserted in the initial configuration at 
p = 0.0765A - 2 , close to the melting density, with the exact zero temperature 
Shadow Path Integral Ground State method. The crystalline order is found to be 
stable also in presence of many vacancies and we observe two completely different 
regimes. For small n„, up to about 6, vacancies form a bound state and cause a 
decrease of the crystalline order. At larger n v , the formation energy of an extra 
vacancy at fixed density decreases by one order of magnitude to about 0.6 K. 
In the equilibrated state it is no more possible to recognize vacancies because 
they mainly transform into quantum dislocations and crystalline order is found 
almost independent on how many vacancies have been inserted in the initial 
configuration. The one-body density matrix in this latter regime shows a non 
decaying large distance tail: dislocations, that in 2D are point defects, turn out 
to be mobile, their number is fluctuating, and they are able to induce exchanges 
of particles across the system mainly triggered by the dislocation cores. These 
results indicate that the notion of incommensurate versus commensurate state 
loses meaning for solid 4 He in 2D, because the number of lattice sites becomes ill 
defined when the system is not commensurate. Crystalline order is found to be 
stable also in 3D in presence of up to 100 vacancies. 
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1. INTRODUCTION 

Solid 4 He is the subject of many experimental and theoretical studies[TJ El [3] since 
the discovery of non classical rotational inertia^ (NCRI), an expected manifestation 
of supersolidity. [5 Supersolidity is a new state of matter in which spatial order and 
off-diagonal long range order (ODLRO) are present at the same time, implying some 
form of superfluid properties, and was already predicted long ago.[6l[7] This novel 
state of matter attracts interest also in the prospect of finding it in cold atoms 
in optical lattices. [5] However the precise nature of solid 4 He at low temperature 
is still elusive. Experimental evidence suggests that defects play an important role 
in NCRI, but which kind of disorder is responsible for the anomalous properties of 
solid 4 He is still not clear. Many different defects have been considered, but none 
of the proposed models has shown to be able to capture all the phenomenology of 
supersolidity. For example the stiffening of 4 He below the "transition" temperature[9] 
suggests dislocations as a candidate, but this possibility has difficulty in explaining 
the NCRI seen in solid 4 Hc in Vycor[J or in Aerogel. (TPl Grain boundaries have been 
considered too,[T] but NCRI has been observed also in single crystals. [TT] 

Defected solid He systems have been studied by means of Quantum Monte 
Carlo (QMC) techniques. [H [Il[Tl[Til[Tl[T71[ia[Tl[5ni[lIl[21[5g Such QMC 
methods are by construction equilibrium methods, so defects have to be stabilized 
by periodic boundary conditions (via a suitable choice of the simulation box) or by 
fixing the degrees of freedom of a number of atoms surrounding the interested defect. 
This constraint on the configurational space usually does not prevent the study of 
the physical properties of the defected system, and the results are generally in a 
good quantitative agreement with experimental observations; see, for example, the 
estimation of the vacancy activation energy [12] or the binding energy to the dislocation 
core of a 3 He atom. [35] Exact QMC results, both at finite[T71 HHJ \W\ and zero 
temperature [24], show that 2 and 3 vacancies form a bound state. A systematic study 
of many vacancies in solid 4 He is still lacking; questions like whether the crystalline 
structure is stable in presence of a large number of vacancies, whether the vacancy- 
vacancy interaction is pairwise additive (or almost so), whether many vacancies lead 
to phase separation as suggested in Ref. jTTJ QI5] , or whether vacancies turn themselves 
into other kinds of defects are, according to us, still unanswered. 

Since the earliest days of supersolidity quantum lattice models have been 
considered. [25] A terminology borrowed from lattice models is currently used for solid 
4 He in terms of a commensurate state (i.e. a perfect crystal in which there is an 
integer number of atoms per unit cell) or an incommensurate one (a crystal with non 
integer probability of occupation of the unit cell). QMC at T = K of solid 4 He with 
one or few vacancies in the simulation cell beautifullv|16) confirms the picture of an 
incommensurate state: the vacancy(-ies) is (are) mobile and the number of density 
maxima is equal to the number of particles plus the number of vacancies. But are we 
sure that we understand vacancies at a finite concentration in a macroscopic system? 
Is the option commensurate or incommensurate really appropriate for solid 4 He? One 
should keep in mind that in a lattice model or in cold atoms in an optical lattice 
the lattice periodicity is externally imposed, whereas in solid 4 He the same dynamical 
entities, the atoms, have to build up the periodicity as well as the derealization 
required for supersolidity. Such "mermaid" aspect of the atoms is unique to a quantum 
solid. The purpose of this article is to address such issues with exact QMC methods 
at T = K in two dimensions (2D). 
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We find that size and box commensuration effects are very pronounced, so in 
order to be able to explore very large distances in a Monte Carlo calculation, up to 
100 A, we have mainly studied solid 4 He in 2D. In this context, the 2D solid is also of 
interest as a simple model for out of registry solid 4 He adsorbed on planar substrates, 
such as graphite or silica glasses. New experimental investigations [26] are under way 
to answer the question whether a supersolid state is present also in adsorbed 4 He. We 
have studied 2D solid He systems at T — K with a number of atoms up and above 
thousand containing up to 30 vacancies. We find that crystalline order is stable also 
in presence of a large number of vacancies, at least in the range < 2.5% of vacancy 
concentration studied by us. Multiple vacancies are highly correlated with a preference 
to form linear structures, but in presence of 10 or more vacancies the scars in the lattice 
are healed away transforming the vacancies mainly into dislocations. Such dislocations 
are not permanent structures, but are mobile and their number is fluctuating. This 
allows exchange of particles across the system, not only in the cores of dislocations, 
giving the possibility of establishing a well defined, at least local, phase [17] and indeed, 
in the simulated systems, we find a non decaying large distance tail in the one-body 
density matrix. This hints at the possibility that an extended system is supersolid 
but more work is needed to corroborate this conclusion. 

We do not address here the origin of vacancies/dislocations, i.e. if the true ground 
state of solid 4 He contains or not zero-point defects. Our finding is that solid 4 He, 
at least in 2D, is much more complex than the picture derived from the notion of 
commensurate vs incommensurate state because the number of lattice sites is an 
ill defined quantity when dislocations are present. The number of lattice sites is a 
meaningful quantity, at least for the sizes we are able to simulate, only if the number 
of particles and the simulation box are such that the regular lattice exactly fits in or the 
misfit is limited to just one or very few vacancies. We have performed few simulations 
also in three dimensions (3D). Also in this case we find that, near melting, crystalline 
order is stable even in presence of 100 vacancies, but we have not yet performed a 
detailed characterization of the disorder in the system. 

The paper is organized as follows: Sec. [2] deals with the exact T = K SPIGS 
method. Details on the simulations are outlined in Sec. [3] Sec. 0] contains our results 
and our conclusions are given in Sec. [5] 

2. THE SPIGS METHOD 

We employ the exact T = K Shadow Path Integral Ground State (SPIGS) [TS] 
method, an "exact" QMC technique. [28, 29 SPIGS is an extension of the Path Integral 
Ground State (PIGS) [3D] method. The aim of PIGS is to improve a variationally 
optimized trial wave function ipT by constructing a path in the Hilbert space of the 
system which connects the given ipT to the lowest energy state of the system, ipo, 
constrained by the choice of the number of particle N, the geometry and the boundary 
conditions of the simulation box and the density p. During this "path" , the correct 
correlations among the particles arise through the "imaginary time evolution operator" 
e~ rH , where H is the Hamiltonian operator. For a large enough r, an accurate 
representation for the lowest energy state wave function is given by tp r = e~ TH ipT, 
which can be written analytically by discretizing the path in imaginary time and 
exploiting the factorization property e~( T1+T2 ) if = e~ TlH e~ T2H . In this way, ip T turns 
out to be expressed in terms of convolution integrals which involve the "imaginary time 
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propagator" (R\e~ \R') for small St, for which very accurate approximations can 
be found in literature. [3lJ E2] This maps the quantum system into a classical system of 
open polymers. [30] An appealing feature peculiar to the PIGS method is that, in ip T , 
the variational ansatz acts only as a starting point, while the full path in imaginary 
time is governed by e~ rH , which depends only on the Hamiltonian operator. We 
have recently shown that PIGS results for large enough r are unaffected by the choice 
of ipT both in the liquid and in the solid phase [2"5l 129) thus providing an unbiased 
exact T = K QMC method. Within SPIGS a shadow wave function (SWF) [331 HO 
is taken as tpT and this choice was shown to greatly accelerate convergence to tpQ. 
Another appealing feature of the SPIGS method is that it recovers the solid phase via 
a spontaneously broken translational svmmetrv|15j and there is no constriction on the 
atomic positions, so it is particularly indicated in studying crystals with defects. 

3. MODEL SYSTEM AND SIMULATION DETAILS 

By studying the equation of state we find the freezing and the melting densities 
to be, respectively, pf — 0.0672 h~ 2 and p rn — 0.0724 A~ 2 . As simulation box 
we choose then a rectangular box fitting a regular triangular lattice with M lattice 
sites at p = 0.0765 A~ 2 . Periodic boundary conditions (pbc) are applied in both 
directions. Dealing with low temperature properties, 4 He atoms are described as 
structureless zero-spin bosons, interacting through the HFDHE2 Aziz potential. [55] 
In order to avoid any approximation associated with the estimation of tail corrections 
due to the finite size of the simulation boxes, we have considered a truncated and 
shifted Aziz potential which goes to zero at r cu t = 6 A. Since r cu t is well below the 
minimum size of each considered box, the results relative to different box sizes can be 
compared without any further correction. Our SPIGS computation is based on the 
pair-product approximation for the imaginary time projector [3 1) and the imaginary 
time step St — 1/40 K _1 has been chosen in order to ensure a good accuracy and a 
reasonable computational effort. 28 

An important parameter in the computation is the total projection time r. Any 
trial wave function ipT can be expanded in terms of the exact ground state and of the 
excited states of zero total linear momentum. If we call Ae m i n the minimum excitation 
energy for such excited states, r has to be larger than (Ae m i n ) _1 in order to project 
out any excitation contribution in the rate of rejection of such contributions 
depends also on the quantity one is computing due to sensitivity to missing long 
range correlations in ipx or to wrong short range correlations. In the perfect crystal, 
using K as energy unit, r = 0.775 K _1 was found [2"B] to be large enough to ensure 
convergence of both diagonal and off-diagonal observables when a SWF is used as 
ipT- This is compatible with the criterion r > (Ae m j n ) -1 . In fact, in a perfect 
crystal, the lowest energy excitations are phonons; thus Ae m i n corresponds to the 
combination of two phonons of opposite k with the smallest achievable k value, i.e. in 
a box with pbc fc m i n = 27r/L max where L max is the largest side. Via a new analytic 
continuation method [35] we have estimated the longitudinal sound velocity for 2D 
solid 4 He at p = 0.0765A~ 2 to be c = 36.5 ± 2.5 KA. Then the minimum imaginary 
time requested for convergence is r m j n = 1/Ae m i n ~ 0.25 K _1 in the largest simulation 
box considered here. This value is well below the used one r = 0.775 K _1 . Since the 
transverse sound velocity is typically similar to c, the considered r value is enough to 
ensure also the removal of transverse phonon contribution. The situation might be 
different in presence of vacancies, or other kinds of defects, because novel low energy 
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excitation modes could be present. We have performed several tests in order to be 
confident that also for defected crystals the considered value of r is large enough. The 
most straightforward one is to perform simulations with larger r values and verify that 
the expectation values are compatible within the statistical uncertainty. Such a test 
turns out to be very time consuming and in practice it has been feasible only for a 
number of 4 He atom N < 250, where we have increased the value of r up to 1.16 K _1 
finding good convergence for diagonal properties. An indirect test on the convergence 
is provided by the behavior of suitably chosen observables during the imaginary time 
path. If we define 

{tpT'\W2T-T') 

for < t' < t, when r is large enough we have ?/>2t-t' — ipo and Eq. [1] gives the 
mixed expectation value 0(t') — (ip T > \O\1p0) / {ip T i \tpo) as a function of the partial 
imaginary time r'. If 0(t') has a plateau when r' — > r, we are guaranteed that r 
is large enough for convergence and the plateau represents the "exact" expectation 
value of O. A couple of examples are shown in FigQ]for the potential energy (O = V) 
and for the static structure factor at the lowest wave vectors {O — p_fc min /Ofc min , where 
p~ = -j= Y^j=i ex P (ik ■ fj) is the density fluctuation operator) along the simulation 
box axis. As far as the potential energy V is concerned, excitation modes are rapidly 
projected out and r ~ 0.25 K _1 is enough to achieve convergence both in the perfect 
and in the defected systems, as shown in Fig[T]by the presence pf a plateau of V(t') 
in both cases. Larger values of the imaginary time are needed to reach convergence 
for quantities that depend on long-range correlations, such as the static structure 
factor at small wave vectors (5fc min ) .[37] In the perfect crystal, convergence for Sfc min 
is achieved at about 0.6 K _1 , as shown in FigHJ In the defected crystal (especially 
when many vacancies are present as in the case considered in the figure) the chosen 
r = 0.775 K _1 value is just enough to get convergence. An interesting physical effect 
emerges from FiglTJ the static structure factor is strongly affected by the presence of 
defects especially along the nearest neighbor direction. 

Two independent simulations are needed in order to calculate the formation 
energy of n v vacancies, one with N — M particles (the regular ideal crystal, dubbed 
also perfect crystal), and one with n v less atoms (defected crystal). In this second case 
we say that we have n v = M — N vacancies, even if in general we can talk of vacancies 
only in the starting configuration of the simulation, where n v atoms are removed 
from the perfect lattice. In fact, in SPIGS, like in PIMC, there is no constraint on 
the atomic positions, so that vacancies are free to transform into different kinds of 
defects. In presence of vacancies, after removing n v particles from the ideal starting 
configuration, we rescale the dimensions of the simulation box to reset the system to 
the original density. This is performed to circumvent the need of correcting the energy 
due to density changes caused by the inclusion of vacancies. The formation energy of 
n v vacancies is proportional to the difference [HJ [35] between the energy per particle 
of the defected crystal, e(M — n v ), and of the perfect one, e(M), i.e. 

AE nv = (M - n v )[e(M - n v ) - e(M)]. (2) 

When vacancies turn into dislocations this A.E n ^ has the meaning of defect formation 
energy at fixed density. One can consider also the formation energy A.E n ^ at a fixed 
lattice parameter; this is given by 

AE nv = AE nv - n v (j., (3) 
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Figure 1. Upper panel: Mixed expectation value for the potential energy per 
particle as a function of the projection time r' computed in a 2D 4 He crystals at 
p = 0.0765A™ 2 in the M = 960 box, with n v = (open symbols) and n v = 20 
(filled symbols). Lower panel: Static structure factor computed in the same 
systems for the lowest wave-vectors along the first neighbor direction (circles, 
fcmin = 0.055A -1 ) and orthogonal to the first neighbor direction (diamonds, 
^'min = 0.059A - 1 ). Open and filled symbols have the same meaning ad in upper 
panel. 



fj, being the chemical potential, that for our system turns out to be [i — 14.22 ±0.02K. 

We check the presence of crystalline order by monitoring the static structure factor 
for the presence of Bragg peaks. In addition, we compute the particle coordination 
number, via a Delaunay triangulation [39] (DT) of the sampled configurations, in 
order to estimate the amount of disorder in the system. [3D] In a perfectly ordered 
2D triangular crystal each atom is linked to 6 other atoms in the DT. Atoms with 
number of coordination not equal to 6 are then a measure of the disorder in the 
crystal. [40] A conservation law exists for the coordination numbers: 40 
5 00 

^(6-0^ = ^(^-6)^ (4) 

i=3 i+7 

where Ni is the number of i— coordinated atoms, so that we can consider only with 
i < 6. Therefore we take Xd = 2N4 + N$ as an estimate of disorder in the system (we 
never found 3-sided polygons, i.e. N3 = 0). In a perfect (N = M) 2D quantum crystal 
the coordination is 6 only on the average, in fact, due to the large zero point motion, 
fluctuations are present such that atoms do not have always coordination 6. Thus, a 
more useful quantity measuring the net amount of disorder in a crystal with defects is 
the difference between the observed Xd and that of the corresponding perfect crystal: 

X d - (2N 4 + N 5 )lj - (2N 4 + N 5 )ir°- (5) 
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Figure 2. Defect formation energy A_E n „ as a function of riu in 2D solid 4 He at 
p = 0.0765A" 2 computed in boxes with different lattice site numbers M. Dashed 
lines are linear fit to the data. 



When the number of vacancies n v is small, up to 6, we have also computed the 
vacancy- vacancy correlation function g vv (r) by recording the relative distances among 
vacancies during the Monte Carlo sampling. The determination of the vector positions 
of the vacancies in a crystalline configuration is far from being trivial due to the large 
zero point motion, to high vacancy mobility and because in our algorithm the center 
of mass is not fixed. Vacancy positions are obtained via the coarse-graining procedure 
illustrated in Ref. .24 . This method is efficient as long as we have few vacancies. For 
large n v the efficiency in recognizing the position of vacancies becomes extremely poor 
(in fact, as we shall see, vacancies turn into dislocations), then is no more possible to 
compute g V v{r). Different definitions of vacancy positions, like the one in Ref. [18 , 
give very similar results. [41"] 

A special care has been devoted in ruling out possible metastability effects. We 
have considered different starting configurations for vacancies, i.e. we have removed a 
number n v of atoms forming a compact cluster or a linear cluster, or the removed atoms 
come from random lattice positions. After long enough equilibration (our simulations 
are never shorter than 6 x 10 5 Monte Carlo steps) we find agreeing results for systems 
with the same n v value. In order to rule out possible finite size effects and pbc bias 
we have considered systems of increasing sizes (M values). Again we find compatible 
results for systems of different sizes and with equal n v value, so no appreciable finite 
size effects have been detected. The only exception is the one-body density matrix 
that will be discussed below. Furthermore we have considered also crystals with no 
principal crystalline axis parallel to the box sides obtaining once more the same results. 
Moreover, we have considered systems where the starting configuration was obtained 
by removing n v particles from an equilibrated configuration of the perfect crystal and 
at the same time the positions of one (two) line(-s) of atoms, parallel to one (both) 
box side(-s) were kept fixed during the Monte Carlo sampling. Even in this cases, the 
results do not change, for instance, for large n v values, vacancies are found to turn 
into dislocations with the same features as in the fully mobile systems. 
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Figure 3. Vacancy-vacancy correlation function gw(r) computed in 2D solid 
4 He at p = 0.0765A -2 in a box with M = 240 lattice sites for different number of 
vacancies n v . g V v(r) is normalized with the correlation function of n v vacancies 
randomly distributed on the same lattice. 
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Figure 4. AS m ^/N = 5Sax/(M - n„) - Smax /M as a function of the number 
of vacancies n v (and of the concentration of vacancies x v in the inset) for different 
M. Sjiax is the main Bragg peak integrated intensity. [43 j The main Bragg peak 
is at \k\ ~ 1.87A - 1 . Lines are aid to the eyes. 



4. RESULTS 

We have studied a triangular crystal at p = 0.0765A" 2 with M = 240, 480, 960 and 
1440 lattice positions of the perfect crystal. Our results for the formation energy 
AE„ v at constant density are plotted in Fig. [5J The dependence of AE n ^ on n v 
is monotonic and systematically sublinear, confirming the existence of an attractive 
interaction among vacancies. Systems with the same n v and different M have the same 
A.E nv within the statistical uncertainty, i.e. AE nv has no significant dependence on 
x v = n v /M, at least in the range x v < 2.5 % that we have studied. From the plot in 
Fig. [5] it is possible to recognize two different behaviors: for n v < 6, AE Uv deviates 
rapidly from the linear dependence and vacancies form a bound state as indicated 
by a (roughly) exponentially decreasing correlation function g vv {r), plotted in Fig. [3] 
For n v > 10 the ratio AE nv /n v remains practically constant with a value of about 
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0.61 K. This means that, when some vacancies are already present in the system, the 
creation of an additional vacancy has a very low cost, about one tenth of the cost of 
a single vacancy. Notice that in the linear regime A.E nv , at fixed lattice parameter, is 
negative, there is a gain of 13.61 K for each additional vacancy. 

The double regime behavior of AE nv is reflected also in other properties of the 
system, for instance in the static structure factor S(k). For a finite crystal the height 
•Smax of the Bragg peaks has a strong dependence on the number of particle N and 
Sm&x/N shows a slow convergence (like 7V-2/3 in 3D02] and N^ 1 / 2 in 2D) to the 
thermodynamic limit. Note that being at T = K, the crystalline order is stable 
also in 2D. Our computation of S^^/N in the perfect crystal verifies this iV -1 / 2 
dependence as expected for a 2D crystal. This TV -1 / 2 dependence of S^^/N arises 
from missing phonon modes for k < fc m i n in a finite box so that essentially the same 
TV -1 / 2 contribution is expected to be present in the defected system. Therefore by 
taking the difference between S max /N in the defected system and the value in the 
prefect crystal such TV -1 / 2 contribution cancels out to a large extent so we can compare 
results for crystal of different sizes and n v . Since the crystalline lattice relaxes around 
a vacancy, one expects that the heights of the Bragg peaks in S(k) decrease with 
increasing number of vacancies. This is observed, as shown in Fig. |4l only when n v 
is small. For large n v values we find a very different behavior. The main Bragg peak 
height has some broadening but its integrated intensity 43 does not show a significant 
dependence on n v , rather it oscillates around an almost constant value. In Fig. [3] we 
show the static structure factor S(k) in the k plane obtained for n v = 0, 5, 15 and 25 in 
the largest considered box with M = 1440 lattice sites. When the number of vacancies 
is small (Fig. [5^) we find Bragg peaks in the first reciprocal vector star as sharp as 
in the case of the perfect crystal, within the resolution Ak = due to the finite 
size of the simulation box, but with a decreased height. Otherwise, see Fig. [SJ: and d, 
for large n v , the Bragg peaks turn out to be slightly broadened, but their integrated 
intensity (see Fig. [4j shows no significant dependence on n v . The broadening cr max of 
the main Bragg peak has a linear dependence on the concentration of defects x v , with 

"max — 1 .5 X X v A 

A similar behavior is displayed also by the amount of disorder Xd, plotted in 
Fig. [U Xd increases with n v for small n v values and saturates to an almost constant 
value for larger n v . By examining sampled configurations of the system in this 
regime, we find that it is no more possible to identify vacancy positions. Rather 
one finds dislocations, as shown by a DT of the position of the atoms, see Fig. [TJi 
and b. Dislocation cores in 2D are point defects characterized by a couple 5-7 in 
the coordinations number of neighboring atoms in the DT. We find that dislocations 
do not prefer any of the principal lattice directions. Studying a large collection of 
independent configurations one finds that a large fraction of them contains a couple 
of dislocation cores (see Fig. [7J. For example, for n v = 20 in the M — 960 box, there 
is a probability of about 90% for two dislocation cores and 10% of finding more than 
two dislocation cores. In addition there is a probability around 3% of finding also 
some isolated vacancies besides dislocations (see Fig. [7]:). Similar values are found for 
n v = 20 in the M = 1440 box. We stress that, even if initially placed in a compact 
cluster configuration, vacancies do not separate into a vacancy rich region surrounded 
by a regular crystal. Rather they mostly reorganize themselves across the system 
giving rise to dislocations. In this regime of large n v the quantity that maintains 
direct physical meaning is N = M — n v , the number of atoms. When atoms rearrange 
themselves giving rise to dislocations, the number of lattice sites becomes ill defined 
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Figure 5. (color online). Static structure factor S(k) in the k x ,k y plane 
computed in a 4 He 2D crystal at p = 0.0765A -2 in a box with M = 1440 lattice 
sites and with different n v . The wave vectors are in A -1 . 




Figure 6. Amount of disorder as a function of n v . Lines are aid to the eyes. 
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Figure 7. (color online). Dclaunay triangulation of typical equilibrium 
configurations of a 2D 4 He crystal at p = 0.0765A - 2 with (a) n v = 10 and 
M = 480, (b) n v = 10 and AI = 490 and (c) n v = 20 and M = 1440. The 
bold(red) lines indicate the interrupted lines of atoms. We report the coordination 
number for atoms only when different from 6. The couple 5-7 indicates a 
dislocation core in 2D and the clumps of two 5 and two 7 (shadowed area in 
panel a) indicate bound pairs of dislocations. [40] In (c), besides the dislocations, 
also two close threefold symmetric monovacancies 44 are present (shadowed area) . 
Notice that in all panels only a small part of the full system is shown. 
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Figure 8. (color online). Vacancy configuration in a 2D 4 He with M = 240 and 
(a) n v = 2, (b) n v = 4 and (c) riu = 6. Vacancy positions are obtained as in 
Ref. 1241 . Notice that in all panels only a small part of the full system is shown. 



because the lattice positions of the perfect crystal compatible with the simulation box 
and the pbc do no more represent equilibrium positions for the particles, and it is 
no more possible to recognize a regular lattice describing the equilibrium positions of 
the particles. Still we continue to use M and n v = M — N as a convenient measure 
of deviation from the ideal crystal N = M where crystalline order is perfect. In 
the regime of small n v , vacancies are strongly correlated and prefer to form linear 
configurations, as shown in Fig. [51 These linear configurations can be considered 
as forerunners of the quantum dislocations found at larger n v . This behavior has 
some similarities with the behavior of colloidal crystals in 2D.[44l [45] We find that 
dislocations are not fixed structures, rather they are very mobile and their number 
changes along the simulation. As an example, in Fig. [H] we show the evolution of the 
DT of a configuration for n v = 10 vacancies in the box with M — 480 sampled at 
different Monte Carlo times. 
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Figure 9. (color online). Example of the evolution of DT at different Monte Carlo 
times (MCS) in a converged computation for a 2D 4 He crystal at p = 0.0765A - 2 
with n v = 10 and M = 480. The bold (red) lines indicate the interrupted lines of 
atoms. We report the coordination number for atoms only when different from 6. 
The couple 5—7 indicates a dislocation core in 2D. The clump of two 5 and two 7 
(shadowed areas) indicates a bound pair of dislocations (BPD). 40 Starting from 
the configuration shown in the panel a, after 300 MCS dislocation cores are found 
to be glided along the horizontal lattice axis by a lattice parameter (b). After 
1100 MCS (c) we find three dislocation cores. Then after 1700 MCS we find again 
two dislocation cores, aligned along the vertical simulation box axis, plus three 
BPD. In all panels the whole system is shown. 

The one-body density matrix pi(r,r') has a central role since the presence of a 
plateau at large \r — r ' indicates the presence of ODLRO, i.e. the system has Bose- 
Einstein condensation (BEC). Even px(r,f) seems to show the signature of a double 
regime behavior depending on the number of vacancies. Some of our results for the 
one-body density matrix are reported in Fig. 1101 As found in the 3D case|16j. when 
only a single vacancy is present pi displays a non-zero plateau, which is the signature 
of ODLRO. As vacancies are added, the tail of one-body density matrix is depressed, 
even if it does not show a simple exponential decay. In order to be conclusive on the 
behavior of p\ in the large distance range one should compute it in larger systems; 
we have faced these computationally intensive calculations in the n v range where 
vacancies give rise to quantum dislocations. Here (sec Fig. ITOb and Fig. Wlk ) p\ 
displays a recognizable plateau at large distances (very evident in the largest system) . 
In order to understand the origin of the observed off-diagonal contributions to p\ in the 
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Figure 10. One-body density matrix pi computed along the first neighbor 
direction in 2D 4 He crystals at p = 0.0765A - 2 in different boxes in (a) vacancy 
regime and (b) dislocation regime. 




Figure 11. (color online). Delaunay triangulation of typical equilibrium 
configurations of a 2D 4 He crystal at p = 0.0765A -2 in an off-diagonal SPIGS 
simulation with n v = 20 and M = 960. The bold (red) lines indicate the 
interrupted lines of atoms. We report the coordination number for atoms only 
when different from 6. The couple 5-7 indicates a dislocation core in 2D. The 
clump of two 5 and two 7 indicates a bound pair of dislocations. 1401 The blue dots 
are the snapshots of the half-polymers. [15] The two half-polymers are mainly 
found to occupy the dislocation cores like in (a), but they have been found also 
outside the core (both half-polymers in (c) or only one in (b)). Notice that in all 
panels only a small part of the full system is shown. 

dislocation regime it is useful to recall that in a SPIGS computation pi(r, r') (which is 
the probability amplitude to destroy a particle in f ' and to create one in r) is obtained 
by splitting one of the linear polymers in two half polymers. |15) one departing from r 
and the other from r '. We find that the main contribution to the plateau comes from 
configurations where at least one of the half-polymers occupies a dislocation core, see 
Fig. UlTa.b'). This means that the system is able to transfer particles from a quantum 
dislocation to another. This is quite distinct from the observed superfluidity along the 
core of a screw dislocation[21] in 3D, in fact in 2D a dislocation core is a point defect 
and not a linear one as in 3D. Furthermore, the possibility of finding a half-polymer 
out of the cores like in Fig. ITTT b.c) means that quantum dislocations are able also to 
induce vacancies in the surrounding crystal. Contrarily to diagonal properties, in the 
computed one-body density matrix we find presence of significant size effects, as can 
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Figure 12. (color online). One-body density matrix pi (r, r') in the x — x'.y — y' 
plane computed in a defected 4 He 2D crystal at p = 0.0765A -2 . 

be inferred for example by the significant difference between the pi in the M — 700, 
n v = 10 system and in the M — 960, n v = 20 one, in the intermediate distance range 
(20-40A). This is shown more clearly, in Fig. [T2lwhere we report p\ in the full plane, 
one can notice the presence of ridges in the tail of pi in the smaller systems (M = 480 
and M — 700, see panels (a) and (b)) that are no more present in the largest system 
(see Fig. IT2"b). Such ridges are due to a commensuration effect with the simulation 
box. We have evidence that in the M — 480 and M = 700 systems dislocation cores 
prefer to stay at a distance of about (n v ±2) x a (a is the lattice parameter) that turns 
out to be comparable with L x /2. This is no more true in the box with M = 960, in 
which case p\ has no ridge in the tail (see Fig. [T2h ) . 

In order to determine how the plateau of p\ is affected by finite size effects one 
should perform computations for an even larger number of particles, bur N ~ 10 3 is the 
maximum number that presently we can handle. Another reason for studying larger 
systems is to establish how the number of dislocations scales with the size of the system 
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and with the number of vacancies. Presumably supersolidity and BEC are present in a 
macroscopic system only if the concentration of dislocations is finite. In addition there 
is the question if the projection time r is large enough to have convergence to the exact 
result. On the basis of indirect tests performed on diagonal observables which depend 
on long-range correlations (see Fig. [IJ we are rather confident that the used value 
of t is large enough to get convergence, but we are not able to provide a direct test 
by using a larger value of r again for computational limitations (a computation with 
r = 0.775K" 1 with St = 1/40K" 1 and 10 3 atoms has to handle about 2 x 10 5 degrees 
of freedom). Our results in the dislocation regime hint that ODLRO is present, but 
further tests are needed in order to firmly establishing this. In any case our results 
on off-diagonal properties in 2D systems enlighten a novel ability of dislocations in 
inducing exchanges even in their surroundings. 

5. CONCLUSIONS 

In summary, we find that crystalline order in 2D 4 He is stable also in presence of a large 
number of vacancies and there is no tendency to phase separation. This is true also 
in 3D, we have studied up to 2548 particles and 98 vacancies and in all cases we find 
that crystalline order is present, as shown by the presence of Bragg peaks in the static 
structure factor. When, in the 2D system, the number of vacancies is of order of 10 
and above, vacancies inserted in the initial configuration loose their identity and most 
of them become quantum dislocations. Such dislocations turn out to be very mobile 
and are able to induce exchange of particles across the system, which is necessary 
to supersolidity. The ability of dislocations to induce vacancies in the surrounding 
crystal could be relevant also for the 3D case. If this feature will be confirmed in 3D, 
the superfluidity would be not restricted to the dislocation cores, [3T] but exchange 
processes necessary for supersolidity would be promoted by vacancies even in the 
bulk crystal far from dislocation cores. Simulations with numbers of particle orders of 
magnitude larger are needed to answer the question whether the number of dislocations 
cores will increase with the system size and with the number of vacancies to give a 
finite concentration of dislocations able to induce ODLRO even in a macroscopic 2D 
4 He crystal, but this is actually beyond our computational possibilities and we have to 
leave it for future investigations. Moreover, how much the phase correlation triggered 
by dislocations depends on the concentration of dislocations and its relevance for the 
3D solid 4 He are still open questions. Our results mainly indicate that the usual 
concepts of commensurate or incommensurate solid, borrowed from lattice models, is 
not appropriate for solid Helium, a system where the crystal lattice is not externally 
imposed, but is self-induced by the correlation among particles. In fact, the number 
of lattice sites M is found to be an ill defined quantity when the crystal houses many 
vacancies/dislocations. 

We do not address the origin (intrinsic or extrinsic) of the vacancies/dislocations 
studied in the present paper. The issue if the ground state could contain zero point 
defects is still debated. [TH [3j |46j [47] In case the ground state of solid 4 He contains 
defects, our findings suggest that the Andreev-Lifshitz-Chester scenario [BJ [7] would in 
a certain sense revive, but in terms of ground state dislocations rather than vacancies. 
On the other hand, even if the ground state has no zero point defects, our findings 
suggest that in presence of extrinsic disorder this will manifest more in terms of 
fluctuating dislocations rather than vacancies, at least at low T and in 2D. 
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